program main
   implicit none
   real::x,y(1),h
   real::xa,xb
   integer::i
   xa=0.0
   xb=2.0
   h=(xb-xa)/10
   x=xa
   y(1)=0.5
   write(*,"(3A13)")"x","exact","rk4"
   do i=1,10
      call rk4(1, h, x, y, derivs)
      write(*,*)x,exp(-x)*(1.5*exp(x/2)-1),y
   end do
contains
   subroutine derivs(x,y,k)
      real,intent(in)::x
      real,intent(in)::y(*)
      real,intent(out)::k(*)
      k(1)=(exp(-x)-y(1))/2
   end subroutine derivs
end program main

subroutine rk4(n, h, x, y, derivs)
   implicit none
   integer,intent(in)::n
   real,intent(in)   ::h
   real,intent(inout)::x,y(n)
   interface
      subroutine derivs(x,y,k)
         real,intent(in)::x
         real,intent(in)::y(*)
         real,intent(out)::k(*)
      end subroutine derivs
   end interface
   real:: hh
   real::k1(n),k2(n),k3(n),k4(n)
   hh = h/2
   call derivs(x    , y       , k1)
   call derivs(x+hh , y+hh*k1 , k2)
   call derivs(x+hh , y+hh*k2 , k3)
   call derivs(x+h  , y+h*k3  , k4)
   x =x+h
   y =y+(k1+2*k2+2*k3+k4)*h/6
end subroutine rk4

